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I. 



INTRODUCTION 



Simulations of Quantum Chromodynamics (QCD) are nowadays being carried out at almost physical parameters. 
Therefore studies of hadron structure within lattice QCD are beginning to yield results that can be connected to 
experiment more reliably than ever before. In these lattice QCD studies one calculates matrix elements of local oper- 
ators between hadron states. Unless these operators correspond to a conserved current they have to be renormalized. 
Calculation of renormalization factors can be carried out using lattice perturbation theory. Although perturbation 
theory on the lattice is computationally more complex than in the continuum these calculations can be extended 
beyond one-loop order [1-3] . Various methods to improve the convergence of lattice perturbation theory have been 
introduced [4, 5] yielding valuable first input to the values of the renormalization constants. In this work we will 
use the perturbative results to improve the non-perturbative evaluation of the renormalization constants. We use 
the Rome-Southampton method (also known as the RI-MOM scheme) [6] to compute renormalization coefficients of 
arbitrary quark- antiquark operators non-perturbatively. In this approach the procedure is similar to that used in 
continuum perturbation theory. In particular, the renormalization conditions are defined similarly in perturbative and 
non-perturbative calculations. The renormalization factors, obtained for different values of the renormalization scale, 
are evolved perturbatively to a reference scale /x = 2 GeV. In addition, they are translated to MS at 2 GeV using 3-loop 
perturbative results for the conversion factors. Since in the end one wants to make contact with phenomenological 
studies, which almost exclusively refer to operators renormalized in the MS scheme of dimensional regularization, one 
needs the renormalization factors leading from the bare operators on the lattice to the MS operators in the continuum. 

A number of lattice groups are producing results on nucleon form factors and first moments of structure functions 
closer to the physical regime both in terms of pion mass as well as in terms of the continuum limit [7-13]. In 
these lattice QCD computations one calculates hadron matrix elements of bilocal operators. In order to compare 
hadron matrix elements of local operators to experiment one needs to renormalize them. The aim of this paper 
is to calculate non-perturbatively the renormalization factors of the vector, axial-vector, scalar, pseudoscalar and 
tensor currents within the twisted mass formulation of Wilson lattice QCD [14]. We show that, although the lattice 
spacings considered in this work are smaller than 1 fm, 0(a 2 p 2 ) terms are non-negligible and are significantly larger 
than statistical errors. We therefore compute the 0(a 2 g 2 )-terms perturbatively and subtract them from the non- 
perturbative results. This subtraction suppresses lattice artifacts considerably depending on the operator under study 
and leads to a more accurate determination of the renormalization constants. This approach was applied to evaluate 
the renormalization constants for one-derivative bilinear operators in Ref. [15]. 

The paper is organized as follows: in Section II we give the expressions for the fermion and gluon actions we 
employed, and define the operators. Sections III and IV concentrate on the perturbative procedure, and the 0(a 2 )- 
corrected expressions for the renormalization constants Z q and Zq. In Section V we provide the renormalization 
prescription of the RI'-MOM scheme, and we discuss alternative ways for its application, while in Section VI we 
provide all necessary formulae for the conversion to MS and the evolution to a reference scale of 2 GeV. Section VII 
focuses on the non-perturbative computation, where we explain the different steps of the calculation. The main results 
of this work are presented in Section VIII: the reader can find numerical values for the Z-factors of the fermion field 
and fermion operators, which arc computed non-perturbatively and corrected using the perturbative 0(a 2 ) terms 
presented in Sections III and IV. For comparison with phenomenological and experimental results, we convert the 
Z-factors to the MS scheme at 2 GeV. In Section IX we give our conclusions. 



Our perturbative calculation makes use of the twisted mass fermion action including the usual clover (SW) term 
with a clover parameter that is left free. For Np flavor species and using standard notation, this action reads 



II. FORMULATION 



A. 



Lattice actions 




■+a pi>f (x + an)+ ip f (x + a n) (r + 7^ ) U x+a ^ x 



ipf(x) 



x,f,H 





(1) 



x, f, n, v 
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where the Wilson parameter r is henceforth set to r = 1, / is a flavor index, = [7^, 7^/2 and the clover coefficient 

csw as well as the twisted mass parameter Mo are kept as free parameters throughout. is the standard clover 
discretization of the gluon field tensor [16]. 



For the non-perturbative calculation, we consider the purely twisted mass fermion action (no clover term), which 
for two degenerate flavors of quarks is given by: 

Sf = a 4 X(*)(^(^ + %) - Y^^l +m °+ W75T 3 )x(x) , (2) 

X 

with t 3 the Pauli matrix acting in the isospin space, and Mo the bare twisted mass. Maximally twisted Wilson quarks 
are obtained by setting the untwisted bare quark mass mo to its critical value m cr , while the twisted quark mass 
parameter mo is kept non- vanishing in order to give the light quarks their mass [14, 17]. The physical quantities 
extracted from lattice simulations employing maximally twisted Wilson quarks are automatically 0{a) improved [17]. 
In Eq. (2) the quark fields \ are in the so-called "twisted basis" . The "physical basis" is obtained for maximal twist 
by the simple transformation: 

ip{x) = exp (y75T 3 ^ x(x), i>(x) = x(x) exp (y75T 3 ^ , w = | . (3) 

In terms of the physical fields the action is given by: 



St = a 4 ]T V(*) + - ^75T 3 (-| + m cr ) 



Mo UW- W 



One can check that this action is equivalent to the action in the twisted basis given by Eq. (2), by performing the 
rotations defined in Eq. (3) and identifying mo = m cr . 



For gluons we employ the Symanzik improved action, involving Wilson loops with 4 and 6 links (lxl plaquette, 
1x2 rectangle, 1x2 chair, and lxlxl parallelogram wrapped around an elementary 3-d cube), which is given by 



Sg = ~o 

9o 



plaq. 



rect. 



chair paral. 



(5) 



The coefficients Cj can in principle be chosen arbitrarily, subject to the following normalization condition, which 
ensures the correct classical continuum limit of the action: 



c + 8ci + 16c 2 + 8c 3 = 1. (6) 

Some popular choices of values for Cj used in numerical simulations will be considered in this work, and are itemized 
in Table I (the acronym TILW represent the Tadpole Improved Liischer-Weisz action); they are normally tuned in 
a way as to ensure 0(a 2 ) improvement in the pure gluon sector. In our non-perturbative computation presented 
here we employ the tree- level Symanzik action (co = 5/3, c\ = —1/12, c 2 = C3 =0). Our 1-loop Feynman diagrams 
do not involve pure gluon vertices, and the gluon propagator depends only on three combinations of the Symanzik 
parameters: 

Co = c + 8ci + 16c 2 + 8c 3 = 1, 

Ci=c 2 +c 3 , (7) 
C 2 = ci - c 2 - c 3 . 



Therefore, with no loss of generality we can set c 2 = 0. 
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Action 






Cl 


C3 


Plaquette 




1.0 








Symanzik 




5/3 


-1/12 





TILW, /3c = 


8.60 


2.3168064 


-0.151791 


-0.0128098 


TILW Br-n — 


8.45 


2 3460240 


-0 154846 


-0 01 34070 


TILW, /3c Q = 


8.30 


2.3869776 


-0.159128 


-0.0142442 


TILW, Pc = 


8.20 


2.4127840 


-0.161827 


-0.0147710 


TILW, /3c = 


8.10 


2.4465400 


-0.165353 


-0.0154645 


TILW, /3c = 


8.00 


2.4891712 


-0.169805 


-0.0163414 


Iwasaki 




3.648 


-0.331 





DBW2 




12.2688 


-1.4086 






TABLE I: Input parameters Co, ci, C3. 



B. Definition of operators and Renormalization condition 

The ultra-local bi-fermion operators considered in this work are the following: 



0% 



= Xr a X 



= X75i ' X = 



O a v = 



= Xlv T X 



?pT a ip a = 1,2 
— i'0751'0 a = 3 

'^ 75 T a V> a = 1,2 
—iipltp a = 3 



= X7s7^°X =< 



^lbltiT 2 ip a = 1 

- - < - - 0757^ T V a = 2 
^j^t 3 ^ a = 3 

V>7^r 2 V> a = 1 
-Vh^t^ a = 2 
k V'757m 1-3 V 7 a = 3 
^nuT a tp a = 1,2 
-#75(7^^1^ a = 3 

^^^T a i\) a = 1,2 
—i'tpa^ttp a = 3 



(8) 
(9) 

(10) 

(11) 

(12) 
(13) 



For convenience we have included Of* even though its components are related to those of 0f>. We denote the 
corresponding Z-factors by Zg, Zp, Z v , Z A , ZJ, Z^ p . In a massless renormalization scheme the renormalization 
constants are defined in the chiral limit, where iso-spin symmetry is recovered. Hence Z-factors become independent 
of the isospin index a = 1,2,3 and we drop the a index on the Z-factors from here on. Still note that, for instance, 
the physical V^r 1 ^ is renormalized with Za while f/>7 M T 3 ^> needs Zy, which differ from each other even in the chiral 
limit. 



The renormalization constants are computed both perturbatively and non-perturbatively in the RI'-MOM scheme 
at different renormalization scales. We translate them to the MS scheme at fi =2 GeV using a conversion factor 
computed in perturbation theory to 0(g 6 ) as described in Section VI. The Z-factors are determined by imposing the 
following conditions in the massless theory, i.e., at critical mass and vanishing twisted mass 



(14) 



^l2 Tr 



= 1. 



(15) 
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where the trace is taken over spin and color indices, /j, is the renormalization scale, while S L and T L correspond to 
the perturbative or non-perturbative results, and is the tree-level result for the propagator defined as: 

s<°> M = -wy . (16) 

while r^°) is the tree-level value for the fermion operators S, P, V, A, T, T', that is 

r (0) (f>) = l, 75, 7/*, 75 7m> 75<V/, a^, (17) 

respectively. The trace is taken over spin and color indices. For alternative renormalization prescriptions the reader 
can refer to Ref. [15]. 

The choices for and given in Eqs. (16) - (17) arc optimal, since we obtain Z q = 1, Zq = 1 when the 
gauge field is set to unity. Similarly, in the perturbative computation this condition leads to Z q = 1, and Zq = 1 at 
tree-level. 

III. CORRECTIONS TO THE FERMION PROPAGATOR 

The fermion propagator of the interacting theory is given by the following 2-point correlation function (Green's 
function), with the various quantities computed in perturbation theory: 



(x% f (x)x b p 9 (y)) = J " j^s ab e^-y^s tree ■ ]T (-s^ p (p) • s tree ) n j , (is) 



with S^^p) being the amputated, 1PI 2-point function in momentum space, computed perturbatively up to a desired 
order. Stree is the tree-level propagator for the twisted mass action and it is given by 

M(p) = ml + - ^sin 2 (ajJ M /2), $ = ^ 7/J - sin(ap M ) , (19) 



Stree — , ■ , , 

if + M (p) + ^ 7 5 r 3 a ~ ~ a 



where a, /3 are Dirac indices, /i, fi are flavor indices in the fundamental representation of SU(Np), and a, b are color 
indices in the fundamental representation of SU(N C ). The dot product runs over flavor and Dirac indices. Due to the 
diagonal form of the t 3 matrix, and since we are studying the case of only two degenerate quarks (up/down) we can 
simplify the expression of Stree and omit r 3 by giving a flavor index to the twisted mass parameter, and at the same 
time we take — > m n : 

Stree = ~ 7J\ > (20) 

ip + M(p) + ifi ( J V 

where now fi^ = +/i for the up quark propagator and ^ = — for the down quark propagator. The 1-loop 
Feynman diagrams that enter our 2-point amputated Green's function calculation (S^ ), are depicted in Fig. 1. 





1 2 

FIG. 1: One-loop diagrams contributing to the fermion propagator. 
Wavy (solid) lines represent gluons (fermions). 
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For the algebraic operations involved in evaluating Feynman diagrams, we make use of our symbolic package in 
Mathematica. In a nutshell, the required steps for the computation of a Feynman diagram are the following (the 
reader can find more details in Ref. [18]): 

• A preliminary expression for each diagram can be obtained by contracting the appropriate vertices, which is 
performed automatically once the algebraic expression of the vertices and the topology ( "incidence matrix" ) of the 
diagram are specified. To limit the proliferation of the algebraic expressions we exploit symmetries of the theory, and 
we simplify the color dependence, Dirac matrices and tensor structures. 

• The 0(a 2 ) computation introduces several complications, especially when isolating logarithms and Lorentz non- 
invariant terms, which leads to a whole family of infrared divergent integrals. These can be reduced to a minimal set 
of approximately 250 'basis' integrals. This is achieved by converting all propagator denominators to a standard form 
(q 2 + M 2 )^ 1 using two kinds of subtractions, one for the fermion propagator 

1 1 r4E AJ sin 4 (^/2)-4(E M sin 2 (^/2)) 2 -4moE^in 2 (< ZM /2) 
q 2 q 2 | q 2 q 2 

where the denominator of the fermion propagator, q 2 , is defined as 

( Z 2 =^sin 2 (^)+(mo + ^ 2 ) 2 + ^ /) ) 2 , 9 2 = 4^sin 2 (|), M 2 = m 2 + M 2 , (22) 

and one for the gluon propagator : 

D(q) = D plaq (q) + {D(q) - D plaq (q)} = D plaq (q) + D plaq (q) {^(g) - D" 1 ^}^) . (23) 

D is the 4x4 Symanzik gluon propagator; the expression for the matrix (^D p ia q (q) — D^ 1 (q)j, which is 0(q i ), is 
independent of the gauge parameter, A, and it can be easily obtained in closed form. Moreover, we have 

(D plaq ( q )), v = 6 f-(1-X) q ^ (24) 

Terms in curly brackets of Eqs. (21) and (23) are less IR divergent than their unsubtracted counterparts, by one 
or two powers in a. These subtractions are performed iteratively until all divergent integrals (initially depending 
on the fermion and the Symanzik propagator) are expressed in terms of the gluon propagator, (q 2 + M 2 ) -1 . The 
computation of the divergent integrals is performed in a non-integer number of dimensions D > 4. Ultraviolet 
divergences are explicitly isolated a la Zimmermann and evaluated as in the continuum. The remainders are D- 
dimensional, parameter-free, zero external momentum lattice integrals which can be recast in terms of Bessel functions, 
and finally expressed as sums of a pole part plus numerical constants. 

A small subset of the infrared divergent integrals, shown in Appendix A, contains the most demanding ones in the 
list; they must be evaluated to two further orders in a, beyond the order at which an IR divergence initially sets in. 
As a consequence, their evaluation requires going to D > 6 dimensions. A correct way to evaluate strong divergent 
integrals is given in detail in a previous publication [18]. 

• The required numerical integrations of the algebraic expressions for the loop integrands are performed by highly 
optimized Fortran programs; these are generated by our Mathematica 'integrator' routine. Each integral is expressed 
as a sum over the discrete Brillouin zone of finite lattices, with varying size L (4 4 < L 4 < 128 4 ), and evaluated for all 
values of the Symanzik coefficients listed in Table I. 

• The last part of the evaluation is the extrapolation of the numerical results to infinite lattice size. This procedure 
entails a systematic error, which is reliably estimated, using a complex inference technique; for one-loop quantities we 
expect a relative error smaller than 1CT 7 . 

Next, we provide the total expression for the inverse fermion propagator S~ elt (p), computed up to 1-loop in 
perturbation theory. Here we should point out that for dimensional reasons, there is a global prefactor 1/a multiplying 
our expressions for the inverse propagator, and thus, the 0(a 2 ) correction is achieved by considering all terms up 
to 0(a 3 ). The most general expression for the inverse propagator appears in the Mathematica file Zfactors.m (see 
Appendix D for notation). In the main text we provide a more compact expression, for special values of the various 
parameters, that is tree-level Symanzik improved gluon action, c$w — 0, Landau gauge (A=0), but we keep the 
Lagrangian mass and the twisted mass parameter general. 
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„_i . k . , ap 2 a 2 i 

S pcvt . = m + t(J,j +i p+ — - 



(i 



+g 2 { -13.0232725(2) ip / +m[ 0.5834586(2) - 3 ln[a 2 A/ 2 + a 2 p 2 } - 



(25) 



3A/ 2 lnfl + £ 



M 2 



i\i-f 8.7100834(2) - 31n[a 2 AP + a 2 p 2 ] 



3A/ 2 ln[l 



M 2 



- 0.8530378(3) m 2 - 1.842911859(4) M 2 
+ i mp 1 (0.3393996(2) 



M 2 



6M 2 m 2 ln[l J- 
, ^hi[a 2 M 2 + a 2 p 2 } 



40.69642965(5) p 2 



3M 2 . 3, , 2 „ /r2 , ^ 3M 4 ln[l + ^] 



2p 2 2 



2(p 2 



>= / , , r o , o „ 6A/ 2 ln[l + Ay] ' 
+ i ^ m 7° -6.68582372(4) + 3 ln[a 2 i\/ 2 + a 2 p 2 \ + l - MIL 



+ a* 



to 2.3547298(l)p 2 + 2.3562747(1) m 2 + 3.46524146(4) M 1 



M 4 M 6 3m 2 p 2 



3m z 



L1M ^ \n[a 2 M 2 + a 2 p 2 } + ( V - - 9m 2 - 2M 2 - 



6p 2 3(p 2 ) 2 2(M 2 +p 2 ) 
M e \ M 2 ln[l + \ 



3(p 



2 ^2 



/ M 4 A/ 6 3to 2 d 2 

;^7 5 0.70640552(8) p 2 + 6.79538844(2) to 2 + 1. 16985307(3) Af 2 -—r 4—^- y „. 

y 6p z 3(jH) 2 2(Az z +fr) 



, 2M ^i„r„2» r2 , „2„ 2l , /V n ^2 ^ 2 ^ 6 ^Ml + l^] 



„2 OA/f2 

V 3m 

4 3 



(^-9m 2 - — - 3 



EpPpJ (m + i/ry 5 ) / 1 2M 2 M 4 2M e / j 2M 6 \Af 2 ln[l 



2 9p 2 3(p 2 ) 2 3(p 2 ) 3 V 3 3(p 2 ) 



p 

TP 



ip 1 1.1471643(7) p 2 - 0.2145514(2) to 2 + 1.15904388(6) M 2 



9to 2 Af 2 209M 4 M 6 7M 8 



2p 2 360p 2 240(p 2 ) 2 40(p 2 ) 3 



/ 73p 2 3to 2 2M 2 



V 360 2 3 



+ — - ln[a 2 Af 2 + aV] + — + — -=- + 



1 9m 2 43M 2 M 4 7M 6 \ Af 4 ln[l + *,] 



24 ' 2p 2 ' 72p 2 12(p 2 ) 2 40(p 2 ) 3 



67M 2 Af 4 8A/ 6 7A/ 8 157, r2 „ r2 , 2l 

^ ^.2478764(2) - ^ + ^-^ - ^ + ^ - ^ ln[a 2 A/ 2 + aV] 

1 5M 2 5M 4 7Af 6 \ M 4 ln[l + ^]\ * (E p ' J^) j/ / 7 Af 2 67A/ 4 



2 18p 2 12(p 2 ) 2 30(p 2 ) 3 y (p 2 ) 2 



+ 



p 2 I 240 48p 2 72 (p 2 ) 2 



13M 6 7Af 8 / 5 5A/ 2 M 4 7M 6 \Af 4 ln[l + ^] 



24(p 2 ) 3 12(p 2 ) 4 V 12 V 4(p 2 ) 2 12(p 2 ) 3 y (p 2 ) 



+ 0(a 3 ,.g 4 ) 
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where $ = X) P 7 P Pp an d "S^ = Sp7 P ? ) p- To make the above expressions less complicated we defined m = mo and 
M 2 = m.Q + Hq. We would like to point out that the up quark propagator is obtained by the choice /i*- 1 -* — +/io, while 
for the down quark propagator one should choose /j,^ = — fi . Moreover, g 2 = \§Ji and Cf = ^£w^ ■ 

Another byproduct of this part of the computation is the additive critical fermion mass; its general expression 
depends on csw and the Symanzik parameters. These are terms proportional to 1/a that have been left out of 
Eq. (25) for conciseness: 

+ -0(.9 4 ) • (26) 
a 

The quantities £m (listed in Table II) are numerical coefficients depending on the Symanzik parameters. 



a 



4n ] + e$ csw + et ] clw 



Action 

Plaquette 
Symanzik 
TILW (8.60) 
TILW (8.45) 
TILW (8.30) 
TILW (8.20) 
TILW (8.10) 
TILW (8.00) 
Iwasaki 
DBW2 



-51.4347118(2) 

-40.44324019(7) 

-34.17747288(3) 

-33.9488671(1) 

-33.6344391(1) 

-33.43979350(6) 

-33.1892274(1) 

-32.87904072(9) 

-26.07292275(7) 

-11.5127475(2) 



13.73313097(5; 
11.94821988(5; 
10.76516514(3 
10.71947605(3; 
10.65632621(4 
10.61705314(7 
10.56629305(3; 
10.50313393(3 

9.01533524(3; 

4.9953066(1) 



5.71513853(1) 

4.662672112(4) 

3.998348778(3) 

3.97345187(1) 

3.939135834(8) 

3.917851255(1) 

3.890401337(1) 

3.856345868(2) 

3.1061330684(3) 

1.351772367(3) 



TABLE II: Numerical results for the coefficients Sm\ £m\ Em (Eq. (26)) for different actions. The systematic errors in 
parentheses come from the extrapolation over finite lattice size, L — > oo. 



IV. CORRECTIONS TO FERMION BILINEAR OPERATORS 

In the context of this work we also study the perturbative 0(a 2 ) corrections to Green's functions of local fermion 
operators (currents) that have the form: 

Or = £ £ (*) r «" P" 4' (*)) • (27) 

z a" 0" 

We restrict ourselves to forward matrix elements (i.e. 2-point Green's functions, zero momentum operator insertions). 
The symbol T corresponds to the following set of products of the Euclidean Dirac matrices: 

r e {s,p,y, a,t,t'} e {i, 7 5 , 7(I , 7 5 7 „ 7 5 v , v } ; ^ = ^,7.], ( 2 ») 

for the scalar 0$, pseudoscalar Op, vector Oy, axial Oa, tensor Ot and tensor prime Ot> operator, respectively. The 
matrix elements of Ot 1 can be related to those of Ot; this is a nontrivial check for our calculational procedure [18]. 
The relationship between the amputated 2-point Green functions At and Kt> is: 

= — — £p v v i A j,/ . (29) 

/J,' v' 

The matrix elements of the above set of fermion bilinear operators can be obtained as: 

(x% f (x)OrX b 9 (y)) = J\ ^5 ab e^-y^S ■ Ap-(p) • S^" (30) 



where Ap Crt '(p) is the amputated 1PI 2-point Green's function of each operator Or, in momentum space, which upon 
contraction of indices becomes: 

An?) = £<*£(P) (xUp)T a >p,x%(p)) X$( P )) amp -- (31) 

a 1 0' 

The only 1-particle irreducible Feynman diagram that enters the calculation of the above Green's function is shown 
in Fig. 2. 




FIG. 2: One-loop diagram contributing to the bilinear operators. A 
wavy (solid) line represents gluons (fermions). A cross denotes the 
Dirac matrix T. 

In this diagram there are two fermion propagators, for which we allowed different ^-values, in order to have more 
general results. In other words, each of the two internal fermion lines on the left and on the right of the operator 
insertion (see Fig. 2) can independently represent the up or down propagator. For the evaluation of the Z-factors, we 
keep the two flavors independent. The amputated Greens functions of the operators Or are given in the Mathematica 
file Zfactors.m. As mentioned above, one may choose the two fermion propagators of the diagram to correspond either 
to the up or down quark, thus there are two twisted mass parameters and ^ 2 \ These can have any sign and 
the only restriction is: \ = Im^'I- 

V. QUARK FIELD AND QUARK BILINEAR RENORMALIZATION CONSTANTS IN THE RI'-MOM 

SCHEME 

An operator renormalization constant (RG) can be thought of as the link between its matrix element, regularized on 
the lattice, and its renormalizcd continuum counterpart. The RCs of lattice operators are necessary ingredients in the 
prediction of physical probability amplitudes from lattice matrix elements. In this section we present the multiplicative 
RCs, in the RI'-MOM scheme, of the quark field (ZP ert ) and quark bilinear operators (Zp Crt ), obtained by using the 
perturbative expressions of S~ 1 (p) and Ap(p). 

The RI'-MOM renormalization scheme consists in imposing that the forward amputated Green function Ar(p), 
computed in the chiral limit and at a given (large Euclidean) scale p 2 — [i 2 , is equal to its tree- level value. In practice, 
the renormalization condition is implemented by requiring that in the chiral limit 1 : 

Z-'Zr V r (p)\ Pp ^ p = 1, Vr( P ) = \tt[At(p) ■ P r ] , (33) 
where Pr are the Dirac projectors defined as follows: 

Pr € {P s ,Pp,Pv,Pa,Pt,Pt>} = {1, l\ -7V -7 5 ^ ~^}\ (34) 

they are chosen to obey the relation Tr[r • Pr] = 4. The traces are always taken over the spin indices. The quark field 
RC Z q , which enters Eq. (33), is obtained by imposing, again in the chiral limit, the condition 2 : 



1 A simpler version of Eq. (33) is given by the relation: 

Z- 1 Z r -Tr\A T (p) • Ap ee l = iTrU^ 00 • AjH , 



(32) 



where Ap ee is the tree- level value of Ap(p). 
2 Strictly speaking, the renormalization condition of Eq. (36) defines the so called RI' scheme. In the original RI-MOM scheme the quark 
field renormalization condition reads: 

7m" 



Z- 1 — Tr 
q 16 



1 . 



(35) 



The two schemes differ in the Landau gauge at the N 2 LO. 
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Z- l V q (p)\ Pp ^ p = l, V q ( P ) 



iEp7pSin(ap P ) 



s~Hp) 



(36) 



We compute Z q in the RI'-MOM renormalization scheme, dehned in Eq. (14), which can be Taylor expanded up 
to 0(a 2 ) terms. This leads to: 



EnJpiPp- \P\ 



a 2 EpPp ' 

3 Z P A 



^l-loop(P) 



■0(aV 



J p p =Pp 



-Tr 



2/ 



1 j/p4 



2 p 2 (p 2 ) 2 



Si-loop (*>) 



(37) 



The trace is taken only over spin indices and S^J^p is the inverse fermion propagator that we computed up to 1-loop 
and up to 0(a 2 ). We make the following definitions for convenience: p 2 = J2 P P% P^ = ^1 P P% = ^plpPp anc ^ 
/ 3 = Ep JpP% 

A very important issue is that the 0(a 2 ) terms in Z q depend not only on \p\, but also on the direction of the 
renormalization scale, p p , as manifested by the presence of J2 P P P : 



vr rt b) 



-Tr 



T 2 EpPp 



Spert.O) 



+ 0(a i g 2 ,g i ). 



(38) 



As a consequence, alternative renormalization prescriptions, involving different directions of the renormalization scale 
fJ-p =Ppi treat lattice artifacts differently. 

By implementing the perturbative expressions of S^ 1 (p) and Ap(p) in Eqs. (33) and (36), we obtain the corre- 
sponding RCs. For the following special choices (independent): tree- level Symanzik gauge action, csw = 0, Landau 
gauge, and general mass m the results of the RCs under study are (for Z q we have also kept fi and M = y/ m 2 + fi 2 
as free parameters): 



report. 1 

g ~ 1 



+g 2 { -13.0232725(2) + am 



0.3393996(2) 



3 ln[a 2 M 2 + a 2 P 2] 3M 2 3M 4 ln[l + 



0.2145514(2) m 2 + 1.15904388(6) M 



2 2p2 2p2 2 

2.1064977(2) pi 9m 2 M 2 209M 4 M 6 



(39) 

1.1471643(7) P 2 



(7ip2 t 3m 2 t 2M 2 t 157p4\ 

180p2 J 



p2 

ln[a 2 M 2 + a 2 p2] 



2p2 



7AP 



360p2 240p2 2 40p2 : 



+ 



M 6 



7AF 



V 360 2 3 

p4 / 43A/ 2 169M 4 
p2 y 80p2 + 180p2 2 + 120p2 3 ~ 20p2 4 

+ 0(a 3 g 2 ,g i ) 



1_ 9m 2 43M 2 M 4 7AI 6 \ Af 4 ln[l + ^] 



24 2p2 72p2 12p2 2 40p2 3 / 



P 2 



1 35A/ 2 Af 4 7A/ 6 \ M 4 ln[l + ^] ' 
12 ~ 36p2 + 6p2 2 + 20p2 3 J p2 2 



Pp =A»p 



The RCs of the bilinear operators have lengthy expressions and are not shown in the main text; they are presented in 
Appendix B. We also include Appendix C which is related to the perturbative results appearing in our publication 
for RCs of one-derivative operators [15]. 
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VI. CONVERSION TO THE CONTINUUM MS SCHEME AT A REFERENCE RENORMALIZATION 

SCALE 

A. Conversion factors 

In this section we provide the expressions for the conversion factors to the MS scheme, as adapted from Ref. [19]. In 
our analysis we use 2-loop formulae; 3-loop corrections for the particular expressions are at the one per cent level. We 
use the following definitions for the conversion factors: Z^ s = C q Z^ 1 _MOM (note we use C q in contrast to Ref. [19] 
where the same quantity was denoted with C" 1 ), and = C X ^'" MOM 



C q = 1 + A^-[(8A 2 + 5)C F 



- (9A 2 - 24C(3)A + 52A - 24((3) + 82) N c + UN F - A 2 C|] % ( 

a \lD7T 



2 



(40) 



C s .p = 1 - (A + 4)C F y|^ + [(24A 2 + 96A- 288C(3) + 57)C F 

+166iV F - (18A 2 + 84A - 432((3) + 1285) N c ] ( ^L- j (41) 

C Ay = l + 0(.g 8 ) (42) 
a 2 

C T t> = l + \C F -jr^ + [(216A 2 + 4320C(3)-4815)Cf- 626ATf 

1671^ LV ' 

+ (162A 2 + 756A - 3024C(3) + 5987) N c ] ^ (yf^) (43) 

The variables g, A correspond to the RI' scheme coupling constant and covariant gauge parameter (defined in Ref. [19]); 
in the Landau gauge, A = 0. ((n) is the Riemann zeta function. The coupling constant, g is related to the bare 
coupling, g , and up to 0(g 6 ) the relation takes the form 

C = |+*(«")(|) 2 +««")(|) 3 - m 

The coefficients d\ and d 2 depend on the renormalization scale a/x and are given by [20] : 

di(a/i) = - ^- ( ^-N c - §JV> ) ln(a/i) - + 2.135730077V C - 0.08414443(8)JV F , 
Zir \ 6 o J ziv c 

ln(aii) 



d 2 {a[i) = d 1 (a^i) 2 



24tt 2 



34\; 2 \, [13N C - A 



3tt 2 

+ t - 2.8626216 + 1.24911587V 2 + N F 

8N? 



0.18898(22) _ 0.15789(26)iV c 



B. Evolution to a reference scale 



All our Z-factors have been evaluated for a range of renormalization scales. In this subsection we use 3-loop 
perturbative expressions to extrapolate to a scale fi = 2 GeV. Thus, each result is extrapolated to 2 GeV, maintaining 
information of the initial renormalization scale at which it was computed. 

The scale dependence is predicted by the renormalization group equation (at fixed bare parameters), that is [21] 

Z$ s M = R fa,H ) )Z% S {l«>) (45) 
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where 



Ro{^^) = 7-^\ (46) 

with 

n*) - £ M») + ^ "f 70 M(A) + ft* + fe* 2 )) 

+ Mto-n-he™-™™ arctan f + \ . (47) 

2/3 pV4ftft ~ ft 2 {y/Wofo - ft J 

To 3 loops, the running coupling, /3- function and anomalous dimension 7 are as follows [21-25], for 7V C = 3: 

gV) = 1 ft ln(ln(^ 2 /A 2 )) 

16^ 2 ftln( M 2 /A 2 ) /3 3 InV/A 2 ) 

+ ^5^3^/^) ln 2 (ln(^ 2 /A 2 )) - ft 2 mQn^/A 2 )) + /3 2 /3 - ft 2 ) + • • • (48) 

Ad = H-|jVf (49) 

OO 

ft = 102 — — Np (50) 
2857 5033 Np 325 

^ = ~2 ~ + ^ (51) 

7o" = (52) 

* = - 2 (t-5"') ( 53 ) 

„ / 20729 79 A/0 . 550 Ar 20 Ar2 \ 

7o S ' P = -2^ (55) 

7f'--2^-^ (56) 



72 



3 9 

2498 + + f C(3))^ + ^f (57) 

lV = lV = lV = (58) 

7 T ' T ' = § (59) 

7i T ' T ' = ~ (26 N F - 543) (60) 

7 2 r,T ' = -J- (36 7V| + 1440 C(3)N F + 5240 iV F + 2784 C(3) - 52555) (61) 
81 
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Eqs. (52) - (61) differ by numerical factors compared to Refs. [21-25] due to alternative definitions of the factor 
#o(a»>Mo)- 

VII. NON-PERTURBATIVE CALCULATION 

In the literature there are two main approaches that have been employed for the non-perturbative evaluation of the 
renormalization constants. They both start by considering that the operators can all be written in the form 

O(z) =J2u(z)J(z,z')d(z'), (62) 

z' 

where u and d denote quark fields in the physical basis and J denotes the operator we are interested in, e.g. JT(z, z') = 
S z ,z'1fi would correspond to the local vector current. For each operator we define a bare vertex function given by 

„12 

G ^) = — E e-^-yHu(x)u(z)J(z,z')d(z')d(y)}, (63) 

x,y,z,z' 

where p is a momentum allowed by the boundary conditions, V is the lattice volume, and the gauge average, denoted 
by the brackets, is performed over gauge- fixed configurations. We have suppressed the Dirac and color indices of G(p). 
The first approach relies on translation invariance to shift the coordinates of the correlators in Eq. (63) to position 
z = [26-28]. Having shifted to z — 0, one can calculate the amputated vertex function for a given operator J for 
any momentum with one inversion per quark flavor. 

In this work we explore the second approach, introduced in Ref. [21], which uses directly Eq. (63) without employing 
translation invariance. One must now use a source that is momentum dependent but can couple to any operator. 
For twisted mass fermions, we use the symmetry S u (x, y) — 75S '(y, x)j$ between the u— and d— quark propagators. 
Therefore with a single inversion one can extract the vertex function for a single momentum. The advantage of this 
approach is a high statistical accuracy and the evaluation of the vertex for any operator including extended operators 
at no significant additional computational cost. Since we are interested in a number of operators with their associated 
renormalization constants we use the second approach. We fix to Landau gauge using a stochastic over-relaxation 
algorithm [29] , converging to a gauge transformation which minimizes the functional 

F = ]T Re tr [U,(x) /})] . (64) 

Questions related to the Gribov ambiguity will not be addressed in this work. The propagator in momentum space, 
in the physical basis, is defined by 

s u (p) e ~ ipix ~ v) (< x Mv)) > sd (p) = yE e ~ wix ~ v) ( d ( x Mv)) • (65) 

x,y x.y 

An amputated vertex function is given by 

T(p) = (S u (p))-'G(p)(S d (p))-'. (66) 
and the corresponding renormalized quantities are assigned the values 

S R (p) = Z q S( P ) , T R {p) = Z- 1 Z V(p) . (67) 

In the twisted basis at maximal twist, Eq. (63) takes the form 

G W = y E e" ip(x_!/) <(1 + h 5 )u(x)u(z)(l + h 5 )J(z, z')(t - i l5 )d{z') d(y)(l - . (68) 

x,y,z,z' 

After integration over the fermion fields, and using S u {x,z) = ^S d \z,x)^$ this becomes 

G(P) = E - h5)S d \z,p)(l - h 5 )J(z, z')(t - i l5 )S d (z',p)(l - i 75 )\ , (69) 
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where (...) G denotes the integration over gluon fields, and S(z,p) = ^ y e zpv S (z , y) is the Fourier transformed prop- 
agator on one of its argument on a particular gauge background. It can be obtained by inversion using the Fourier 
source 

b a a (x) = j* x & aP 8 ab , (70) 
for all Dirac a and color a indices. The propagators in the physical basis given in Eq. (65) can be obtained from 

= j£e-^(i-n5)S d (^)(i-n5)> G 

z 

S » = -\Y, e+iP "(( 1 - i T*)S d \*>P)(*--i'n)) G , (71) 

z 

which evidently only need 12 inversions despite the occurrence of both u and d quarks in the original expression. 

We evaluate Eq. (68) and Eq. (71) for each momentum separately employing Fourier sources over a range of a 2 p 2 
for which perturbative results can be trusted and finite a corrections are reasonably small. 

VIII. NON-PERTURBATIVE RESULTS 

We perform the non-perturbative calculation of rcnormalization constants for three values of the lattice spacing, 
corresponding to = 3.9, 4.05 and 4.20 [30]. In this work we use the lattice spacing as determined from the 
nucleon mass. The values we obtained are 0.089(1) (5) fm, 0.070(1) (4) fm and 0.056(2) (3) fm for j3 = 3.9, 4.05 
and 4.2, respectively [31] and they are in agreement with the ones determined from the pion sector. To extract 
the renormalization constants reliably one needs to consider momenta in the range A^qcd < P < 1/a. We relax 
the upper bound to be ~ 2/ a to 5/ a, which is justified by the linear dependence of our results on a 2 . Therefore, 
we consider momenta spanning the range 0.5 < a 2 p 2 < 5 for which perturbation theory is trustworthy and lattice 
artifacts are still small enough. It is important to note that the extrapolation to the continuum limit, (ap) 2 — > 0, 
is performed for a fixed momentum range in physical units. In Table III we summarize the various parameters 
of the action, that we used in our simulations, and in Table IV we present the values we used for the momenta 
{pt,Px,Py,Pz) = {2n / L t n t ,2-K / L x n x ,2ir / L y n y ,2n/L z n z )). 





a (fm) 


ana 


m v (GeV) 


L 3 


x T 


3.9 


0.089 


0.0040 


0.3021(14) 


24 3 


x 48 


3.9 


0.089 


0.0064 0.37553(80) 24 3 


x 48 


3.9 


0.089 


0.0085 


0.4302(11) 


24 3 


x 48 


3.9 


0.089 


0.01 


0.4675(12) 


24 3 


x 48 


4.05 


0.070 


0.003 


0.2925(18) 


32 3 


x 64 


4.05 


0.070 


0.006 


0.4082(31) 


24 3 


x 48 


4.05 


0.070 


0.006 


0.404(2) 


32 3 


x 64 


4.05 


0.070 


0.008 


0.465(1) 


32 3 


x 64 


4.20 


0.056 


0.002 


0.2622(11) 


24 3 


x 48 


4.20 


0.056 


0.0065 


0.476(2) 


32 3 


x 64 



TABLE III: /^-values and lattice size used in the simulations are given in the first and last columns respectively. The lattice 
spacing a in fm is determined from the nucleon mass. We also give the bare light quark mass a/xo and pion mass. 

The number of configuration in each ensemble varies between 10 to 100. Using even 10 configurations leads to results 
with very high statistical accuracy, easily below 0.5%. Thus, in the plots presented here the statistical errors are too 
small to be visible. We mostly use in our computation democratic momenta, in the sense that they have the same 
p x , p y , p z . We have also tested a few non-democratic momentum, which turn out to behave similarly to democratic 
ones (e.g. (n t , n x , n y , n z )=(3,3,3,2) is similar to (3,3,3,3)). We would like to point out that the non-perturbative 
results have a significant dependence on the value of the momentum in the spatial direction, indicating large lattice 
artifacts in some cases. Such a study appeared in Ref. [15]. 
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P = 3.9 

(n t ,2,2,2), n t : 4-8, 10- 14 
(n t ,3,3,3), th : 2 - 6, 8-9 
(n t ,4,4,4), n t : 4-9 
(3,3,3,2) 



P = 4.05 

(n t ,2,2,2), n t :4-8, 10, 13-14 
(n t ,3,3,3), n t : 2-6, 8-11, 13 
(n t ,4,4,4), n t : 8- 10 
(3,3,3,2) 



g = 4.20 

(n* ,2,2,2), n t :4-8, 10, 13-14 
(n t ,3,3,3), nt : 2 - 6, 8 - 11, 13 
(n t ,4,4,4), n t : 7-11 
(n t ,5,5,5), n t : 2 - 4 
(3,3,3,2) 



TABLE IV: Values of momentum used for the various ensembles at /3 = 3.9, 4.05, 4.20. 



A. Pion mass dependence 



In Table III we give the number of pion mass that we studied for each of the three ft values. These ensembles have 
been produced by the ETM Collaboration [30, 32-34]. The pion mass dependence of Z q , Zy, Za, Zt is displayed 
in Fig. 3 and is not significant. A linear extrapolation to the data shown in Fig. 3 yields a slope consistent with 
zero. This behavior is also observed at the other /3 values. Thus, it would be sufficient to obtain the results on the 
aforementioned RCs at one pion mass value, although we perform the chiral extrapolation with all available data on 
different pion masses. The need of having simulations at a number of pion masses comes from the fact that one has 
to perform the subtraction of the pion pole contribution. This is discussed in Subsection VIII C. 
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FIG. 3: Z q , Zy, Za, Zt at f3 = 3.9, as a function of the pion mass. Computations 
were performed at pion masses of m n = 0.302 GeV (a/io = 0.004), = 0.375 GeV 
(afxo = 0.0064), = 0.429 GeV (afx = 0.0085) and m v = 0.468 GeV (afi = 0.01). 



B. Volume dependence 

We perform the evaluation of the RCs at p = 4.05 and /i = 0.006 for two volumes, 24 3 x 48 and 32 3 x 64 in order 
to check for finite volume effects. For this comparison we used momenta that correspond to the same renormalization 
scale. For the small lattice we use (ap) = 27r(3/48, 3/24, 3/24, 3/24), in lattice units, whereas for the larger one we 
employ (ap) = 27r(4/64, 4/32, 4/32, 4/32). The volume effects appear to be in the worst case ~ 0.1%, as can be seen 
from Table V. We would like to point out that the Zp estimator shows the largest volume dependence, which however 
tends to decrease after the pion pole subtraction (Subsection VIII C, Fig. 4). 



C. Pion-pole subtraction 



The correlation functions of the pseudoscalar operator have pion-pole contamination and therefore need to be treated 
carefully. In order to subtract the pole contribution we use the following Ansatz for the pseudoscalar amputated vertex 
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lattice 




Zs 


Zv 


Zv 


Zk 




24 3 x48 
32 3 x64 


0.82315(7) 
0.82303(3) 


0.743(2) 
0.744(1) 


0.512(2) 
0.508(1) 


0.7068(1) 
0.7069(1) 


0.7935(2) 
0.7935(1) 


0.7759(1) 
0.7759(1) 



TABLE V: Renormalization constants at f5 = 4.05, Ho = 0.006 using two lattice sizes, namely 32 3 x 64 and momentum (4,4,4,4) 
and 24 3 x48 with momentum (3,3,3,3), and at a scale of (ap) 2 ~ 2. 



function, Ap, 

K P = a P + bpml + ^ I , (72) 

which we apply to data produced at a given value of (3. Once we have the fitting parameters we subtract the pion-pole 
using the value of cp, determined from the fitting, i.e. we take 

A s P ub =A P -¥L. (73) 

To reliably obtain the three fitting parameters of Eq. (72) we need the RC of the pseudoscalar operator for at least 
4 pion masses; this is feasible for (3 = 3.9. On the contrary, for (3 = 4.05 we have data for three pion masses, and for 
(3 = 4.20 only for two pion masses. At f3 = 3.9 we determine the parameters using results at three of the four pion 
masses j3 = 3.9 and compare them with the fit resulting when using all available data. The conclusion is that the 
values obtained are compatible. Therefore at (3 = 4.05 we determine the parameters using results at the three pion 
masses that are available. One may observe the effectiveness of the pion-pole subtraction in Fig. 4, where we show 
results before and after the pion pole subtraction. After the subtraction results obtained at different pion mass fall on 
each other. While the pion-pole term has an appreciable contribution, the quadratic term with the bp coefficient is 
expected to be small. The values extracted for b P at (3 = 3.9 and (3 = 4.05 are indeed small showing a very weak pion 
mass dependence of the bp coefficient. This is consistent with the weak pion mass dependence observed for the vector 
and axial- vector RCs (see Subsection VIII A for the other RCs). It is also verified by our data: after subtracting the 
pion-pole term determined from fitting to the data, the remaining pion mass dependence (bp nip) is negligible for all 
the ensembles. This allows us to perform a two parameter fit at (3 = 4.2 of the form: 

Ap = a P + , (74) 
ml 

using data on the two pion masses m v = 476 and 262 MeV. The two sets correspond to different lattice size, 24 3 x48 
(32 3 x 64) for the lower (higher) pion mass. As a result, momenta with the same values for (n t ,n x ,n y , n z ) correspond 
to different (ap) 2 . Thus, in order to perform the fit using the Ansatz of Eq. (74) we carefully choose the momenta 
in the two ensembles to have almost the same (ap) 2 . In general, this could lead to additional uncertainties, but we 
have already checked that volume effects are negligible. Indeed the fit using the Ansatz of Eq. (74) yields a value for 
cp(s) that accurately removes the pion pole as demonstrated in Fig. 5. 



The errors shown in Figs. 4-5 are computed in two ways: using super jackknife error analysis [35, 36] and requiring 
that a correlated change of the fit parameters increases the minimum value of \ by one. We find that both methods 
lead to similar errors. 

Our data for the ratio Zp/Zg also show dependence on the pion mass, and thus to form this ratio we used the 
subtracted data of Fig. 4 which we compute in the chiral limit. This procedure leads to the values shown in Fig. 6. 
With black circles we show the non-perturbative results after subtracting the pion pole from Zp using Eq. (73). 
If one further subtracts from Zs and Zp the perturbative 0(a 2 ) contributions, presented in Sections III - IV, one 
obtains the values shown with the magenta diamonds in Fig. 6. The ratio Zp/Zs is renormalization scale independent 
and therefore one can directly use the C(a 2 )-perturbatively subtracted non-perturbative results to extrapolate to the 
continuum limit, eliminating any remaining cut-off effects. 
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FIG. 4: Zp at /3 = 3.9 (left panel) and /? = 4.05 (right panel) for various masses. The upper plot shows the results before the 
pion pole subtraction as described by Eq. (72), while the lower figure the results upon the appropriate subtraction given in 
Eq. (73). 




FIG. 5: Zp for /3 = 4.20 for the two pion masses. Results are shown upon the pion pole subtraction as described in Eq. (73). 



D. Results in RI'-MOM scheme 



In this section we present results in the RI'-MOM scheme for Z q , Z$, Zp, as well as for the scale- independent RCs 
Zy and Za- We have also performed a computation of Zt and its results are presented in the next section. In all cases 
we subtract the leading discretization effects of 0{a 2 ) computed to one-loop in perturbation theory from the non- 
perturbative results. All renormalization constants are evaluated at the three /3-values, where the simulations were 
carried out. For all /3 values we perform a chiral extrapolation using results at different pion masses; the results have 
negligible dependence on the quark mass as demonstrated in Fig. 3 and therefore we use a constant fit to extrapolate 
to the chiral limit. 

The renormalization constant of the fermion field is needed as an input in various expressions, and the results 
obtained are displayed in Fig. 7. 



The non-perturbative values of Z q are obtained for all available momenta and they reveal a non-smooth behavior 
as a function of the momentum (see Fig. 7), which becomes smoother once we subtract the 0(a 2 ) perturbative 
terms. We would like to point out that in all our non-perturbative results before subtraction we have some data that 
fluctuate several standard deviations around the mean value and these correspond to momenta that lead to large 
non-Lorentz invariant contributions in the perturbative expressions of Sections III-IV. These terms are of the form 
Ep-Pp) / (J2 P Pp)- After subtraction these non-Lorentz invariant contributions are removed (to 0(a 2 g 2 )), resulting in 
the much smoother behavior of the subtracted data. The unsubtracted data of Z q that suffer from large non-Lorentz 
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FIG. 6: Zp/Zs at /3 = 3.9 (upper plot), /3 = 4.05 (middle plot) and /3 = 4.20 (lower plot) as a function of (ap) 2 . In each plot 
we demonstrate the effect of subtracting the C(a 2 )-terms by plotting the non-perturbative results before (black circles) and 
after (magenta diamonds) subtraction. The pion-pole term has been removed from all the data that we show here. 



invariant contributions are show in Fig. 7 with filled black circles. Note that as the lattice spacing decreases, the 
discrepancy between unsubtracted and subtracted data becomes smaller. 

The RCs Zy and Za are scale-independent and therefore there is no need to evolve them. In Fig. 8 we show 
results on Zy and Za before and after subtraction of the perturbatively determined 0(a 2 )-terms. As can be seen, 
the subtraction weakens the dependence on (ap) 2 . In fact, fitting the subtracted data to a straight line of the form 
z + s(ap) 2 results in a value of the slope s consistent with zero. This shows that leading cut-off effects are effectively 
removed by the subtraction of perturbatively determined 0(a 2 )-terms. The small remaining lattice artifacts are 
removed by extrapolating linearly to the continuum line. The unsubtracted data can also be extrapolated linearly 
but, in this case, the slope is generally non-zero as can be seen in Fig. 8. As the lattice spacing decreases, the deviation 
between subtracted and unsubtracted data decreases. We note that, although the value found at the continuum limit 
for the unsubtracted data approaches that extracted for the subtracted data, small differences still remain. This 
is an indication that the systematic error due to cut-off effects is larger than the statistical error and therefore the 
subtraction of C(a 2 )-terms helps in diminishing the uncertainty in the choice of the fit range. 

In order to perform the continuum extrapolation we choose the same momentum range in physical units for all /3 
values and we thus extract all renormalization constants using the same physical momentum range, p 2 ~ 15 — 32 
(GeV) 2 . This momentum range is in line with what has been chosen in our previous work on the RCs for one-derivative 
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FIG. 7: Non-perturbative results on Z q for /3 = 3.9 (upper left plot), /3 = 4.05 (upper right plot) 
and (3 = 4.20 (lower plot). In all plots we show chirally extrapolated results. Black circles (magenta 
diamonds) represent the non-perturbative data before (after) subtracting the C(a 2 )-terms. The 
unsubtracted data that suffer from large non-Lorentz invariant contributions are show with filled 
black circles. 



bilinear operators [15], ensuring that we use data in a region where an approximate plateau exists. The momentum 
range in lattice units at each /3-value is as follows: (3 — 3.9: (ap) 2 = 3 — 5, (3 — 4.05: {ap) 2 = 1.8 — 3.1, (3 = 4.20: 
(ap) 2 = 1.2 — 2.5. The choice for the momentum range is not so relevant for Zy and Za as it was for the case of the 
RCs for one-derivative operators, since the subtracted data are almost constant. However, for consistency we use the 
same range as in Ref. [15]. 



E. MS scheme 



In this Section we present our results on Z q , Z$, Zp and Z^ converted to the continuum MS scheme and at a 
reference scale of fi — 2 GeV. For the conversion from RI'-MOM to MS we use the formulae given in Eqs. (40) - (43). 
We use the 3-loop formulae of Eqs. (45) - (46) to evolve the scale from /x to 2 GeV. 

As already discussed, a "renormalization window" should exist for A.q CD << fi 2 << 1/a 2 where perturbation 
theory holds and finite a artifacts are small, leading to scale- independent results (plateau). In practice such a 
condition is hard to satisfy. The right inequality is extended to (2 — 5)/a 2 leading to lattice artifacts in our results 
that are of 0(a 2 p 2 ). Fortunately our perturbative calculations allow us to subtract the leading perturbative 0(a 2 ) 
lattice artifacts which alleviates the problem. To remove the remaining 0(a 2 p 2 ) artifacts we extrapolate linearly to 
(ap) 2 = as demonstrated in Figs. 8 -11. The statistical errors are negligible, however an estimate of the systematic 
errors is important. The largest systematic error comes from the choice of the momentum range to use for the 
extrapolation to (ap) = 0. One way to estimate this systematic error is to vary the momentum range where we 
perform the fit. Another approach is to fix a range and then eliminate a given momentum in the fit range and refit. 
The spread of the results about the mean gives an estimate of the systematic error. In the final results we give as 
systematic error the largest of the two, which is the one obtained by modifying the fit range. As already mentioned we 
choose the same momentum range in physical units for the three /3-values and extract all renormalization constants 
using the same physical momentum range, p 2 ~ 15 — 32 (GeV) 2 ; within this range the data fall on a straight line of 
a small slope. We note that the 0(a 2 ) perturbative terms which we subtract, tend to decrease with increasing f3, as 
expected. The error bars in Fig. 10 are due to the fit uncertainties in performing the pion pole subtraction. 
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FIG. 8: The renormalization constants for the vector and axial-vector operators, Zv and z?a, for 
/3 = 3.9 (upper left plot), f3 = 4.05 (upper right plot) and /3 = 4.20 (lower plot). In all plots we show 
chirally extrapolated results. Black circles (magenta diamonds) represent the non-perturbative data 
before after subtracting the C?(a 2 )-terms. 
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represent the non-perturbative data before (after) subtracting the 0(a 2 )-terms.The corresponding 
dashed lines show the extrapolation to the continuum limit and the filled diamond shows the final 
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Our final results for the ^-factors in the MS-scheme at \i — 2 GeV are given in Table VI. As pointed out, these are 
obtained in the continuum limit by extrapolating linearly in (ap) 2 using data in a fixed momentum range p 2 « 15 — 32 
(GeV) 2 . The continuum extrapolation was carried out at the chiral limit. The systematic error due to the continuum 
extrapolation, is estimated from the difference between results using the fit range p 2 w 15 — 32 (GeV) 2 and the range 
p 2 « 17 — 24 (GeV) 2 . The results at /3 = 3.9 and j3 = 4.05 agree within error bars with the results of Ref. [26]. Since 
their evaluation procedure differs, this agreement provides a nice confirmation of the values obtained. In Ref. [26] 
the vertex computation employs translation invariance to evaluate the correlation functions for all values of the 
momentum, whereas we calculate the vertex for a given momentum dependent source, leading to smaller statistical 
errors. More importantly, the two procedures differ in the analysis of the lattice data, both in the way the chiral 
extrapolation of the renormalization constants at fixed p 2 is carried out as well as in the way the systematic errors 
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associated with the extrapolation p 2 — > are estimated. Due to the different approaches used, the statistical and 
systematic errors between the two computations is somewhat different. In this work, we additionally compute the 
renormalization constants at {3 = 4.2 but not at (3 — 3.8, which were included in Ref. [26]. Given the consistency 
between our values of the renormalization constant and those of Ref. [26] at (3 = 3.9 and (3 = 4.05 consolidates our 
value at ^ = 4.20 3 . 



p 




Zs 


Z P 


Zp/Zs 


Zv 


Z A 




3.90 
4.05 
4.20 


0.754(9)(9) 
0.775(4)(5) 
0.798(4)(9) 


0.726(5)(11) 
0.691(9)(16) 
0.695(10)(13) 


0.457(10)(16) 

0.497(8)(15) 

0.501(8)(10) 


0.639(3)(1) 
0.682(2)(1) 
0.713(2)(2) 


0.627(1)(3) 
0.662(1)(3) 
0.686(1)(1) 


0.758(1)(1) 
0.773(1)(1) 
0.789(1)(2) 


0.750(9)(10) 

0.798(7)(8) 

0.822(4)(6) 



TABLE VI: Final results of the renormalization constants Z q , Zs, Zp, Zt in the MS scheme, as well as for the scale-independent 
Zp/Zs, Zv and Za- Statistical errors are shown in the first parenthesis. The number in the second parenthesis is the systematic 
error due to the continuum extrapolation, taken as the difference between results using the fit range p w 15 - 32 (GeV) 2 and 
the range p 2 « 1 7 - 24 (GeV) 2 . 



IX. CONCLUSIONS 

The values of the renormalization factor for the fermion field Z q , and for the scalar, pseudoscalar, vector, axial- vector 
and tensor local operators, Z$, Zp, Zy , Z&_, Z^, have been calculated non-perturbatively. The method of choice is to 
use a momentum dependent source and extract the renormalization factors for all the relevant operators. This leads to 
a very accurate evaluation of these factors using a small ensemble of gauge configurations. The precision of the results 
allows us to reliably investigate the light quark mass dependence. For most of the renormalization constants studied 
in this work we do not find any light quark mass dependence within our small statistical errors. For all (3 values we 
obtain the value at the chiral limit by fitting the data to a constant. For the RC of the pseudoscalar operator, Zp, 
we find a quark mass dependence due to the pion-polc, which we subtract. Once the pole is subtracted, the behavior 
of the data show a weak dependence on the light quark mass and therefore we again compute the value at the chiral 
limit by fitting the pion pole subtracted data to a constant. We also show that, despite using a lattice spacing smaller 
than 0.1 fm, 0(a 2 p 2 ) cut-off effects are visible given the high precision with which the RCs are calculated. Thus 
we perform a pcrturbative subtraction of 0{a 2 g 2 ) terms. This leads to a milder dependence of the renormalization 
constants on (ap) 2 . Residual 0(a 2 p 2 ) effects are removed by extrapolating to zero. In this way we can accurately 
determine the renormalization constants in the RI'-MOM scheme. In order to compare with experiment we convert 
our values to the MS scheme at a scale of 2 GeV. The statistical errors are in general smaller than the systematic 
ones. The latter are estimated by changing the window of values of the momentum used to extrapolate to a 2 p 2 = 0. 
Our final values are given in Table VI. 
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3 We note that the preliminary value of Zp = 0.50(2) in MS at 2 GeV was used in quark mass evaluation [37] and b-physics [38] ETMC 
papers is totally consistent with the result of this analysis. 
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Appendix A: Strong IR divergent integrals 

The integrals, with strong IR divergences (convergent only beyond D > 6), encountered in the present calculation 
are listed below with their results. For completeness we also include the integrals that appeared in our related 
publication [15] for the matrix elements of twist-2 operators. All these integrals can be found in electronic form in 
the Mathematica file Zfactors.m, with the names: IntegralPropagatorl - IntegralPropagator3, IntegralBilinearsl - 
IntegralBilinears6, and IntegralExtendedBilinearsl - IntegralExtendedBilinears2. To avoid heavy notation we define: 

M 2 = (mof + Ji 2 , M] = K) 2 + M 2 , 

P 2 = P 4 = Y,P% 

p p 

q u = 2sin(|), ^ 4 ]Tsm 2 (f), 

p 

where q stands for k or k + ap, while k is the loop momentum and p is the external momentum. No summation over 
the indices Vi is implied. 
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Appendix B: Analytic expressions for RCs of bilinear operators 

In this Appendix we provide the analytic expressions for the RCs of the ultra- local bilinears, as defined in Eq. (33): 
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Appendix C: Analytic expressions for one-derivative operators 



Here we present the Z-factors for the one-derivative vector, axial and tensor operators published separately in 
Ref. [15], defined as follows 
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(C2) 



m-Op} _ — , 
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[-iipa^ D p} tip a = 3 
The above operators are symmetrized over two Lorentz indices and are made traceless 
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The one derivative operators fall into different irreducible representations of the hypercubic group, depending on the 
choice of indices: 
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Thus, ^dvij ^dai will be different from Zdv2, ZdA2, respectively. More details on the one-derivative renormalization 
factors can be found in Ref. [15]. 

We have computed, to 0{a 2 ), the forward matrix elements of these operators for general external indices /i, v (and 
p for the tensor operator), external momentum p, m, g, N c , a, cgw and gauge fixing. Our final results were obtained 
for the 10 sets of Symanzik coefficients given in Table I. 

The amputated Greens functions of the Odv operator appear in the Mathematica file Zfactors.m as below: 



^dv — LDV[Action, csw, beta, g2tilde, a, m] + 0(a 3 ,g 4 ), 
^da = LDA[Action, csw, beta, g2tilde, a, m] + C(a 3 , g 4 ), 
^dt = LDT[Action, csw, beta, g2tilde, a, m] + 0(a 3 , g 4 ) 



In order to define Zq, we have used a renormalization prescription which is most amenable to non-perturbative 
treatment: 
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where L° denotes the amputated 2-point Green's function of the operators up to 1-loop and up to 0(a 2 ). These 
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Z-factors appear in electronic form with the name: 
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Due to very lengthy expressions we only show the results for specific choices of the action parameters, that is 
Landau gauge, tree-level Symanzik gluons, csw = 0, m = 0: 
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(z^n)" 1 " 2 = ZDT1[2, 0, 1, g2tilde, a, 0] + 0(a 4 , g 4 ) (C9) 
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where for ZDV1, ZDV2, ZDAl, ZDA2: 
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Appendix D: Notation in Mathematica file: Zfactors.m 

The full body of our results can be accessed online through the Mathematica file Zfactors.m, which is a Mathematica 
input file. It includes the expressions for the amputated Green's functions of the inverse propagator: 

Sport. — ^ 9 ^propagator [Action, csw, beta, g2tilde, a, m, mu] + 0(a \ g 4 )^j , (Dl) 

from which one can construct the fermion field renormalization constant for any renormalization scheme. This 
expression depends on the variables: 

• action: Selection of improved gauge action as follows, 1 — > Plaquctte, 2 — ► Tree Level Symanzik, 3 — > TILW 
{Pc = 8.60), 4 -> TILW (/3c = 8.45), 5 -)• TILW (/3c = 8.30), 6 -» TILW (/3c = 8.20), 7 TILW 
(13 Co = 8.10), 8 TILW (/3c = 8.00), 9 -> Iwasaki, 10 DBW2 

• csw: clover parameter 

• beta: gauge parameter (Landau/Feynman/Generic correspond to 1/0/beta) 

• g2tilde= ^ 67r f g: coupling constant 

• a: lattice spacing 

• m: Lagrangian mass 

• mu: twisted mass parameter 

The expression for the critical mass is defined in the variable mcritical: m cr = mcritical[Action, csw, g2tilde, aL] + 
-0{g 4 ). The reader may also find the amputated Green's functions relevant to the ultra-local operators: 



Ag er ' = scalar [Action, csw, beta, g2tildc, a, m, mul, mu2] + 0(a 3 , g 4 ), (D2) 

Ap Crt ' = pseudoscalar [Action, csw, beta, g2tilde, a, m, mul, mu2] + 0(a 3 ,g 4 ), (D3) 

^pcrt. _ vector [Action, csw, beta, g2tilde, a, m, mul, mu2] + 0(a 3 , g 4 ), (D4) 

^per . _ axial [Action, csw, beta, g2tildc, a, m, mul, mu2] + 0(a 3 , g 4 ), (D5) 

Ay Cr ' = tensor [Action, csw, beta, g2tilde, a, m, mul, mu2] + C(a 3 , g 4 ), (D6) 

Ay° rt ' = tensorprime [Action, csw, beta, g2tilde, a, m, mul, mu2] + (D(a 3 , g 4 ), (D7) 

as well as the Green's functions of the one-derivative vector, axial and tensor operators: 

A^y - = LDV [Action, csw, beta, g2tildc, a, m] + 0(a 3 , g 4 ), (D8) 

A%£ = LDA [Action, csw, beta, g2tildc, a, m] + 0(a 3 , g 4 ), (D9) 

A^*- = LDT [Action, csw, beta, g2tilde, a, m] + 0(a 3 , g 4 ) (D10) 



We note that Eqs. (D2) - (D7) hold for quarks with the same Lagrangian mass and \x\ = ±/Z2, while Eqs. (D8) - 
(D10) correspond to zero twisted mass parameters, /ii, (12, and non-zero Lagrangian mass. 

The RCs are interesting quantities for other studies and are also provided in the Mathematica file Zfactors.m at the 
RI'-MOM scheme by employing Eq. (36) and Eq. (33), for the fermion field and fermion operator RCs, respectively. 
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m, mul, mu2] 


+ 0{a\g 4 ) 


(D13) 




report. 
Z V 


= zv[2, 0,1,1 


g2tilde, a, 


m, mul, mu2] 


+ 0(a 3 ,g 4 ) 


(D14) 




^pcrt. 
Z A 


= za[2,0,l, 


j2tilde, a, 


m, mul, mu2] 


+ 0{a\g A ) 


(D15) 






= zt[2,0,l,j 


*2tilde, a, 


m, mul, mu2] 


+ 0(a 3 ,g 4 ) 


(D16) 


report. \ 
Z Tp ) 


Vl^V2 


= ztp[2,0,l,j 


g2tilde, a, 


m, mul, mu2] 


+ 0(a 3 ,g 4 ) 


(D17) 



The additional variables are 
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•P2: £?=iP? 
•P4: Eti^ 

For completeness we include in the Mathematica file Zfactors.m the conversion factors to the MS scheme for the 
one-derivative RCs, studied in Rcf. [15]. 

Cdvvdai = CDVl[alphaRIprime,lambdaRIprime,CA,Cf,Nf,Tf] + 0(q R ) (D18) 
C D V2,DA2 = CDV2[alphaRIprime,lambdaRIprimc,CA,Cf,Nf,Tf] + 0(g s ) (D19) 

where 

alphaRIprime = g 2 /(16ir 2 ) 
lambdaRIprimc = A = 1 — beta 
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